home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Mac-Source 1994 July
/
Mac-Source_July_1994.iso
/
C and C++
/
Science⁄Math
/
Scientist's Helper src
/
s.helper.5
/
regseries.c
< prev
next >
Wrap
C/C++ Source or Header
|
1986-04-19
|
16KB
|
631 lines
#include "all.h"
#include "regtabext.h"
histogram()
{
int i, n, iCol, *x, *y;
float vMin, vMax, s, *u, binNumber, nm1;
SToR( command.cmdWord[1], &vMin, FALSE);
SToR( command.cmdWord[2], &vMax, FALSE);
if( vMax<=vMin ) {
ErrMsg("bounds out of order");
}
SToI( command.cmdWord[3], &n );
if( (n<1) || (n>table.header.maxRows) ) {
ErrMsg("bad number of bins");
}
s=(vMax-vMin)/n;
nm1 = (float)(n-1);
/* grab temporary memory */
x = NewPtr( (long)(2*n) );
if( MemError()!=noErr ) {
ErrMsg("not enough memory");
}
for( iCol=2; iCol<=table.header.cols; iCol++ ) { /*loop over columns in table*/
/*lock position of column */
HLock( table.ptr[iCol-1] );
/* zero counters */
y=x;
for( i=1; i<=n; i++ ) {
*y=0;
y++;
}
/* count up hits */
u=*(table.ptr[iCol-1]);
for( i=1; i<=table.header.rows; i++ ) {
binNumber = (*u-vMin)/s;
if( (binNumber>=0.0) && (binNumber<=nm1) ) {
*(x+(int)binNumber)+=1;
}
u++;
}
/* copy to table */
u=*(table.ptr[iCol-1]);
y=x;
for( i=1; i<=n; i++ ) {
*u = (float)(*y);
u++;
y++;
}
/*unlock current column*/
HUnlock( table.ptr[iCol-1] );
} /*next column*/
/*dispose of temporary memory*/
DisposPtr( x );
table.header.start = vMin;
table.header.samp = s;
table.header.rows = n;
table.header.interpolated=TRUE;
CreateCol1();
Header2Vars();
}
taper()
{
int xCol, outputCol, startRow, endRow, i;
float f, x, x1, x2, y, y2;
if( (strcmp(command.cmdWord[1],"ascending")!=0) && (strcmp(command.cmdWord[1],"descending")!=0) ) {
ErrMsg( badModifier );
}
SToI( command.cmdWord[2], &xCol );
if( GoodCol(xCol)!=0 ) {
ErrMsg("bad x column");
}
SToI( command.cmdWord[3], &outputCol );
if( (GoodCol(outputCol)!=0) || (table.header.interpolated&&(outputCol==1)) ) {
ErrMsg("bad output column");
}
SToI( command.cmdWord[4], &startRow );
if( (startRow<1) || (startRow>table.header.maxRows) ) {
ErrMsg("bad starting row");
}
SToI( command.cmdWord[5], &endRow );
if( (endRow<1) || (endRow>table.header.maxRows) ) {
ErrMsg("bad ending row");
}
if( endRow<=startRow ) {
ErrMsg("start row, end row out of order");
}
GetTable( startRow, xCol, &x1 );
GetTable( endRow, xCol, &x2 );
f = 3.1415926 / (x2-x1);
if( strcmp(command.cmdWord[1],"ascending")==0 ) {
for( i=startRow; i<=endRow; i++ ) {
GetTable(i, xCol, &x );
GetTable(i, outputCol, &y );
y2 = y * 0.5 * (cos( (double)(3.1415926+f*(x-x1)) ) + 1.0 );
SetTable(i, outputCol, y2, FALSE );
} /*end for i*/
} /*end if ascending*/
else if( strcmp(command.cmdWord[1],"descending")==0 ) {
for( i=startRow; i<=endRow; i++ ) {
GetTable(i, xCol, &x );
GetTable(i, outputCol, &y );
y2 = y * 0.5 * (cos( (double)(f*(x-x1)) ) + 1.0 );
SetTable(i, outputCol, y2, FALSE );
} /*end for i*/
} /*end if descending*/
}
convolve()
{
int inputCol, opCol, n, outputCol, i, j, k, nm1;
float *x, *u, *v, *w, sum;
if( !table.header.interpolated ) {
ErrMsg("table must be interpolated");
}
SToI( command.cmdWord[1], &inputCol);
if( GoodCol(inputCol)!=0 ) {
ErrMsg("bad input column");
}
SToI( command.cmdWord[2], &opCol);
if( GoodCol(opCol)!=0 ) {
ErrMsg("bad operator column");
}
SToI( command.cmdWord[3], &n );
if( (n<1) || (n>table.header.rows) ) {
ErrMsg("bad operator length");
}
SToI( command.cmdWord[4], &outputCol);
if( (GoodCol(outputCol)!=0) || (outputCol==1) ) {
ErrMsg("bad output column");
}
/* grab temporary memory */
x = NewPtr( (long)(4*table.header.rows) );
if( MemError()!=noErr ) {
ErrMsg("not enough memory");
}
/* lock position of columns */
HLock( table.ptr[inputCol-1] );
HLock( table.ptr[opCol-1] );
HLock( table.ptr[outputCol-1] );
/*do convolution*/
u=x;
v=*(table.ptr[inputCol-1]);
w=*(table.ptr[opCol-1]);
nm1=table.header.rows-1;
for( i=0; i<table.header.rows; i++ ) {
sum = 0.0;
for( j=0; j<n; j++ ) {
k = i-j;
if( (k>=0) && (k<=nm1) ) {
sum += (*(v+k)) * (*(w+j));
} /*end if*/
} /*end for j*/
*u=sum;
u++;
} /*end for i*/
/*copy results to table*/
u=x;
v=*(table.ptr[outputCol-1]);
for( i=0; i<table.header.rows; i++ ) {
*v=(*u) * table.header.samp;
u++;
v++;
} /*end for i*/
/* lock position of columns */
HUnlock( table.ptr[inputCol-1] );
HUnlock( table.ptr[opCol-1] );
HUnlock( table.ptr[outputCol-1] );
/*dispose of temporary memory*/
DisposPtr( x );
}
spect()
{
int ifft, n, iCol, i;
float *x, *u, *v, *w;
if( (strcmp(command.cmdWord[1],"amplitude")!=0) &&
(strcmp(command.cmdWord[1],"power")!=0) &&
(strcmp(command.cmdWord[1],"phase")!=0) ) {
ErrMsg( badModifier );
}
if( !table.header.interpolated ) {
ErrMsg("table must be interpolated");
}
n = table.header.rows;
if( n<=128 ) { /*find nearest power of two*/
ifft=128;
}
else if( n<=256 ) {
ifft=256;
}
else if( n<=512 ) {
ifft=512;
}
else if( n<=1024 ) {
ifft=1024;
}
else if( n<=2048 ) {
ifft=2048;
}
else {
ifft=4096;
}
/* grab temporary memory */
x = NewPtr( (long)(4*(ifft+2)) );
if( MemError()!=noErr ) {
ErrMsg("not enough memory");
}
for( iCol=2; iCol<=table.header.cols; iCol++ ) { /*loop over columns in table*/
/* lock position of column */
HLock( table.ptr[iCol-1] );
/* swap into temporary array */
u=x;
v=*(table.ptr[iCol-1]);
for( i=1; i<=table.header.rows; i++ ) {
*u = *v;
u++;
v++;
}
for( i=table.header.rows+1; i<=ifft; i++ ) { /*zero rest of array*/
*u = 0.0;
u++;
v++;
}
cfftr( x, ifft ); /*fast fourier transform*/
u=x;
for( i=1; i<=ifft+2; i++ ) { /*proper normalization*/
*u = (*u) * table.header.samp;
u++;
}
if( strcmp(command.cmdWord[1],"amplitude")==0 ) {
u=x;
v=x+1;
w=*(table.ptr[iCol-1]);
for( i=1; i<=ifft/2+1; i++ ) {
*w = (float) sqrt( (double)((*u)*(*u)+(*v)*(*v)) );
u=u+2;
v=v+2;
w++;
}
}
else if( strcmp(command.cmdWord[1],"power")==0 ) {
u=x;
v=x+1;
w=*(table.ptr[iCol-1]);
for( i=1; i<=ifft/2+1; i++ ) {
*w = (*u)*(*u)+(*v)*(*v);
u=u+2;
v=v+2;
w++;
}
}
else if( strcmp(command.cmdWord[1],"phase")==0 ) {
u=x;
v=x+1;
w=*(table.ptr[iCol-1]);
for( i=1; i<=ifft/2+1; i++ ) {
*w = (float) atan2( (double)(*v), (double)(*u) );
u=u+2;
v=v+2;
w++;
}
}
/*unlock current column*/
HUnlock( table.ptr[iCol-1] );
} /*next column*/
/*dispose of temporary memory*/
DisposPtr( x );
table.header.start = 0.0;
table.header.samp = 1.0 / (ifft*table.header.samp);
table.header.rows = ifft/2+1;
CreateCol1();
Header2Vars();
}
cfour( data, n, isign ) /*fast fourier transform on complex data*/
int n, isign;
float data[];
{
int ip0, ip1, ip2, ip3, i3rev;
int i1, i2a, i2b, i3;
float sinth, wstpr, wstpi, wr, wi, tempr, tempi, theta;
ip0=2;
ip3=ip0*n;
i3rev=1;
for( i3=1; i3<=ip3; i3+=ip0 ) {
if( i3 < i3rev ) {
tempr = data[i3-1];
tempi = data[i3];
data[i3-1] = data[i3rev-1];
data[i3] = data[i3rev];
data[i3rev-1] = tempr;
data[i3rev] = tempi;
}
ip1 = ip3 / 2;
do {
if( i3rev <= ip1 )
break;
i3rev -= ip1;
ip1 /= 2;
}
while ( ip1 >= ip0 );
i3rev += ip1;
}
ip1 = ip0;
while ( ip1 < ip3 ) {
ip2 = ip1 * 2;
theta = 6.283185/( (float) (isign*ip2/ip0) );
sinth = (float) sin( (double) (theta/2.) );
wstpr = -2.*sinth*sinth;
wstpi = (float) sin( (double) theta );
wr = 1.;
wi = 0.;
for ( i1=1; i1<=ip1; i1+=ip0 ) {
for ( i3=i1; i3<ip3; i3+=ip2 ) {
i2a=i3;
i2b=i2a+ip1;
tempr = wr*data[i2b-1] - wi*data[i2b];
tempi = wr*data[i2b] + wi*data[i2b-1];
data[i2b-1] = data[i2a-1] - tempr;
data[i2b] = data[i2a] - tempi;
data[i2a-1] += tempr;
data[i2a] += tempi;
}
tempr = wr;
wr = wr*wstpr - wi*wstpi + wr;
wi = wi*wstpr + tempr*wstpi + wi;
}
ip1=ip2;
}
return;
}
cfftr( x, n ) /*fast fourier transform on real data*/
int n;
float x[];
{
int nn, is, nm, j, i;
int k1j, k1i, k2j, k2i;
float s, fn, ex, wr, wi, wwr, wrr, wwi, a1, a2, b1, b2;
nn = n/2;
is = 1;
cfour( x, nn, is );
nm = nn/2;
s = x[0];
x[0] += x[1];
x[n] = s - x[1];
x[1] = 0.0 ;
x[n+1] = 0.0;
x[nn+1] = (-x[nn+1]);
fn = (float) n;
ex = 6.2831853 / fn;
j = nn;
wr = 1.0;
wi = 0.0;
wwr = (float) cos( (double) ex );
wwi = (float) (-sin( (double) ex ));
for (i=2; i<=nm; i++) {
wrr = wr*wwr-wi*wwi;
wi = wr*wwi+wi*wwr;
wr = wrr;
k1j = 2*j-1;
k1i = 2*i-1;
k2j = 2*j;
k2i = 2*i;
a1 = 0.5*(x[k1i-1]+x[k1j-1]);
a2 = 0.5*(x[k2i-1]-x[k2j-1]);
b1 = 0.5*(-x[k1i-1]+x[k1j-1]);
b2 = 0.5*(-x[k2i-1]-x[k2j-1]);
s = b1;
b1 = b1*wr+b2*wi;
b2 = b2*wr-s*wi;
x[k1i-1] = a1-b2;
x[k2i-1] = (-a2-b1);
x[k1j-1] = a1+b2;
x[k2j-1] = a2-b1;
j -= 1;
}
return;
}
integrate()
{
int i, xCol, yCol, outputCol;
float x1, x2, y1, y2, sum;
SToI( command.cmdWord[1], &xCol);
if( GoodCol(xCol)!=0 ) {
ErrMsg("bad x column");
}
SToI( command.cmdWord[2], &yCol);
if( GoodCol(yCol)!=0 ) {
ErrMsg("bad y column");
}
SToI( command.cmdWord[3], &outputCol);
if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
ErrMsg("bad output column");
}
sum = 0.0;
GetTable( 1, yCol, &y1 );
GetTable( 1, xCol, &x1 );
SetTable( 1, outputCol, sum, FALSE );
for( i=2; i<=table.header.rows; i++ ) {
GetTable( i, xCol, &x2 );
GetTable( i, yCol, &y2 );
sum = sum + 0.5*(y1+y2)*(x2-x1);
SetTable( i, outputCol, sum, FALSE );
x1=x2;
y1=y2;
}
}
differentiate()
{
int i, xCol, yCol, outputCol;
float sum, x1, x2, y1, y2;
SToI( command.cmdWord[1], &xCol);
if( GoodCol(xCol)!=0 ) {
ErrMsg("bad x column");
}
SToI( command.cmdWord[2], &yCol);
if( GoodCol(yCol)!=0 ) {
ErrMsg("bad y column");
}
SToI( command.cmdWord[3], &outputCol);
if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
ErrMsg("bad output column");
}
GetTable( 1, xCol, &x1 );
GetTable( 1, yCol, &y1 );
for( i=1; i<=table.header.rows-1; i++ ) {
GetTable( i+1, xCol, &x2 );
GetTable( i+1, yCol, &y2 );
SetTable( i, outputCol, (y2-y1)/(x2-x1), FALSE );
x1=x2;
y1=y2;
}
SetTable( table.header.rows, outputCol, infinity, FALSE );
}
summation()
{
int i, inputCol, outputCol;
float sum, x;
SToI( command.cmdWord[1], &inputCol);
if( GoodCol(inputCol)!=0 ) {
ErrMsg("bad input column");
}
SToI( command.cmdWord[2], &outputCol);
if( (GoodCol(outputCol)!=0) || ((outputCol==1)&&table.header.interpolated) ) {
ErrMsg("bad output column");
}
sum = 0.0;
for( i=1; i<=table.header.rows; i++ ) {
GetTable( i, inputCol, &x );
sum = sum + x;
SetTable( i, outputCol, sum, FALSE );
}
}
bandpass()
{
int i, inputCol, outputCol;
float flow, fhigh, digt, *fbout, *t, *v;
if( !table.header.interpolated ) {
ErrMsg("cant bandpass uninterpolated table");
}
digt = 1.0 / table.header.samp;
SToR( command.cmdWord[1], &flow, FALSE);
SToR( command.cmdWord[2], &fhigh, FALSE);
SToI( command.cmdWord[3], &inputCol);
if( (GoodCol(inputCol)!=0) || (inputCol==1) ) {
ErrMsg("bad input column");
}
SToI( command.cmdWord[4], &outputCol);
if( (GoodCol(outputCol)!=0) || (outputCol==1) ) {
ErrMsg("bad output column");
}
/* grab temporary memory */
fbout = NewPtr( (long)(4*table.header.rows) );
if( MemError()!=noErr ) {
ErrMsg("not enough memory");
}
/* bandpass input */
HLock( table.ptr[inputCol-1] );
t=*(table.ptr[inputCol-1]);
FILT( t, fbout, table.header.rows, 0, table.header.rows-1, digt, flow, fhigh, 2 );
HUnlock( table.ptr[inputCol-1] );
/* swap result back into table */
HLock( table.ptr[outputCol-1] );
t=fbout;
v=*(table.ptr[outputCol-1]);
for( i=1; i<=table.header.rows; i++ ) {
*v = *t;
v++;
t++;
}
HUnlock( table.ptr[outputCol-1] );
/*dispose of temporary memory*/
DisposPtr( fbout );
}
FILT( FBIN, FBOUT, NN, K, M, DIGT, FLOW, FHIGH, N)
float FBIN[], FBOUT[], FLOW, FHIGH, DIGT;
int NN, K, M;
{
/* CHEBYSHEV RECURSIVE DIGITAL BANDPASS DIGITAL FILTER */
/* FBIN: INPUT ARRAY TO FILTERED */
/* FBOUT: RESULTANT FILTERED ARRAY */
/* NN: LENGTH OF FBIN, FBOUT */
/* K, M: DATA FILTERED BETWEEN THESE LIMITS */
/* IDIGT: SAMPLING RATE */
/* FLOW: LOW PASS FREQ */
/* FHIGH: HIGH PASS FREQ */
/* N: ORDER OF FILTER (MUST BE EVEN OR IT WILL BE MADE EVEN) */
float D[6], G[6], PI=3.1415926, ES=1.0, E, FL, FH, CF, BW, FAC, AYE, BEE;
float PIN, ARG, SR, SI, TEMP;
int N2, I, J;
E=0.1; /* ten percent ripple */
if( (N&1)!=0 ) N=N+1; /*only even orders allowed*/
FL=2.0*FLOW/DIGT;
FH=2.0*FHIGH/DIGT;
/* CALCULATE BANDWIDTH AND CENTER FREQUENCIES IN FRAME WARPED BY Z-FORM */
CF = 4.0 * tan( (double) (FL*PI/2.0) ) * tan( (double) (FH*PI/2.0) );
BW = 2.0 * ( tan( (double) (FH*PI/2.) ) - tan( (double) (FL*PI/2.) ) );
FAC = pow( sqrt((double)(1.0+1.0/(E*E))) + 1.0/E, (double)(1.0/N) );
AYE = 0.5*(FAC-1.0/FAC);
BEE = 0.5*(FAC+1.0/FAC);
PIN = PI/(2.0*N);
N2=N/2;
for( I=1; I<=N2; I++ ) { /* LOOP THROUGH PAIRS OF CONJUGATE POLES */
ARG = (2*I-1)*PIN+PI/2.0;
SR = AYE * cos( (double)ARG);
SI = BEE * sin( (double)ARG);
ES = sqrt( (double)(SR*SR+SI*SI) );
TEMP= 16.0 - 16.0*BW*SR + 8.0*CF + 4.0*ES*ES*BW*BW - 4.0*BW*CF*SR + CF*CF;
D[1] = 1.0;
D[2] = 4.0*(-16.0 + 8.0*BW*SR - 2.0*BW*CF*SR + CF*CF)/TEMP;
D[3] = (96.0 - 16.0*CF - 8.0*ES*ES*BW*BW + 6.0*CF*CF)/TEMP;
D[4] = (-64.0 - 32.0*BW*SR + 8.0*BW*CF*SR + 4.0*CF*CF)/TEMP;
D[5] = (16.0 + 16.0*BW*SR + 8.0*CF + 4.0*ES*ES*BW*BW + 4.0*BW*CF*SR + CF*CF)/TEMP;
TEMP = 4.0*ES*ES*BW*BW/TEMP;
G[1] = TEMP;
G[2] = 0.0;
G[3] = -2.0*TEMP;
G[4] = 0.0;
G[5] = TEMP;
CONVL(FBIN,FBOUT,NN,K,M,D,G);
} /*end for i*/
}
CONVL(FBIN,FBOUT,NN,K,M,D,G)
float FBIN[], FBOUT[], D[], G[];
int NN, K, M;
{
int I, J, IX;
float TEMP;
for( I=K; I<=M; I++ ) { /*SUM PAST VALUES OF INPUT*/
TEMP=0.0;
for( J=1; J<=5; J++ ) {
IX=I-J+1;
if( IX>=0 ) {
TEMP=TEMP+FBIN[IX]*G[J];
}
} /*end for J*/
for( J=2; J<=5; J++ ) { /* SUM PAST VALUES OF OUTPUT */
IX=I-J+1;
if( IX>=0 ) {
TEMP=TEMP-FBOUT[IX]*D[J];
}
} /*end for J*/
FBOUT[I]=TEMP;
} /*end for I*/
}